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Aerodynamic sound of flow past an airfoil 

By Meng Wang 
1. Motivation and objectives 

The long term objective of this project is to develop a computational method 
for predicting the noise of turbulence-airfoil interactions, particularly at the trailing 
edge. We seek to obtain the energy-containing features of the turbulent boundary 
layers and the near-wake using Navier-Stokes Simulation (LES or DNS), and then to 
calculate the far-field acoustic characteristics by means of acoustic analogy theories, 
using the simulation data as acoustic source functions. 

Two distinct types of noise can be emitted from airfoil trailing edges. The first, a 
tonal or narrowband sound caused by vortex shedding, is normally associated with 
blunt trailing edges, high angles of attack, or laminar flow airfoils. The second source 
is of broadband nature arising from the aeroacoustic scattering of turbulent eddies 
by the trailing edge. Due to its importance to airframe noise, rotor and propeller 
noise, etc., trailing edge noise has been the subject of extensive theoretical (e.g. 
Crighton & Leppington 1971; Howe 1978) as well as experimental investigations 
(e.g. Brooks & Hodgson 1981; Blake & Gershfeld 1988). 

A number of challenges exist concerning acoustic analogy based noise computa- 
tions. These include the elimination of spurious sound caused by vortices crossing 
permeable computational boundaries in the wake, the treatment of noncompact 
source regions, and the accurate description of wave reflection by the solid surface 
and scattering near the edge. In addition, accurate turbulence statistics in the flow 
field axe required for the evaluation of acoustic source functions. 

Major efforts to date have been focused on the first two challenges. To this end, 
a paradigm problem of laminar vortex shedding, generated by a two dimensional, 
uniform stream past a NACA0012 airfoil, is used to address the relevant numerical 
issues. Under the low Mach number approximation, the near-field flow quantities are 
obtained by solving the incompressible Navier-Stokes equations numerically at chord 
Reynolds number of 10 4 . The far-field noise is computed using Curie’s extension to 
the Lighthill analogy (Curie 1955). An effective method for separating the physical 
noise source from spurious boundary contributions is developed. This allows an 
accurate evaluation of the Reynolds stress volume quadrupoles, in addition to the 
more readily computable surface dipoles due to the unsteady lift and drag. The 
effect of noncompact source distribution on the far-field sound is assessed using an 
efficient integration scheme for the Curie integral, with full account of retarded-time 
variations. The numerical results confirm in quantitative terms that the far-field 
sound is dominated by the surface pressure dipoles at low Mach number. The 
techniques developed are applicable to a wide range of flows, including jets and 
mixing layers, where the Reynolds stress quadrupoles play a prominent or even 
dominant role in the overall sound generation. 
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2. Accomplishments 

2.1 Aeroacoustic theory 

The density fluctuations due to acoustic wave propagation from a stationary 
aerodynamic source region into a uniformly moving medium is governed by, in 
dimensionless form (Goldstein 1976), 
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is the Lighthill stress tensor defined in terms of the fluctuating velocity relative to 
the free-stream value, v, = u,- — S tl , and 
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represents the viscous part of the Stokes stress tensor. The velocity, density, and 
pressure are nondimensionalized relative to the undisturbed free-stream quantities 
U'oo, p'ooi an d p'ooU'^, respectively. The spatial coordinates are normalized by the 
airfoil chord (or more generally, the characteristic body size) C ' . The time is nor- 
malized by C' /U'oq. Re and M denote respectively the free-stream Reynolds number 
based on chord and the free-stream Mach number. S tJ is the Kronecker delta, and 
the usual summation convention applies for repeated subscripts. 

Like the Lighthill equation, (1) is an exact restatement of the mass and momen- 
tum conservation equations for a compressible fluid. The use of relative velocity in 
the source function ensures that the Lighthill stress, predominantly the fluctuating 
Reynolds stress, is quadratically small outside the source region in the free-stream. 
One notices that, since the radiated acoustic field has a characteristic spatial scale 
of M~ 1 times the hydrodynamic length scale, the two spatial derivative terms on 
the left-hand side of (1) are of O(M) and 0(1), respectively, relative to the time 
derivative term. 

Ffowcs Williams and Hawkings (1969) derived a general solution for noise pro- 
duced by a rigid surface moving through a quiescent medium. An exact solution to 
(1) is most easily obtained by rewriting the Ffowcs Williams-Hawkings equation in 
terms of the “reception coordinates”, i.e., in a reference frame moving with the body. 
For low Mach number flows, however, the bulk convective effect can be ignored to 
first approximation, and the simpler solution owing to Curie (1955) prevails. If x 
and y are used to denote the position vectors of an observation point and a source 
element, respectively, and let r = x - y and r = |r|, Curie shows that 
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for a rigid body at rest. In the above equation p,j = pSij — Tij, and n< is the 
directional cosine of the outward unit normal to the rigid surface S over which 
the surface integration takes place. The volume integral is taken over the entire 
unsteady flow region V external to the body. In the acoustic far-field defined by 
r l e /M, where l e is the typical eddy size, (4) can be simplified to a form most 
suitable for numerical evaluation, 

M 3 d f r 

p ( X)t ) _ i ^ J^n jPij { y,t - Mr)d?y 

+ ^Lwj v r ^ Tii{y ' t - Mr)dZy - (5) 

Furthermore, if the size of the source region is small compared with one acoustic 
wavelength ( l 3 <C l e /M), the source region can be considered acoustically compact. 
The far-field density can be approximated by 

M 3 T: . Ilf 4 rp . . tl 

p(x,t) - 1 « — |^. Di(t - M |x|) + — -j^j - M|x|), (6) 

where 

Dflt) = §ij s n jPtj{y , f)^ 2 y, (7) 
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2.2 Exit boundary correction 

In the application of Curie’s integral solution, the surface integral, taken over the 
finite airfoil surface, is well defined. The same cannot be said, however, regarding 
the volume integral, as the unsteady flow region is often truncated in the wake by 
the artificial computational boundary (cf. Fig. 1). At the outflow boundary the 
Lighthill stress terms are still significant, and their sudden termination are known 
to cause strong, spurious acoustic sources (Crighton 1993). The same difficulty has 
been encountered by Mankbadi et al. (1994) and Mitchell et al. (1995a) in jet noise 
calculations. Mitchell et al. employed model extensions which allow the Lighthill 
source terms to decay to zero slowly downstream of the computational domain. In 
the subsequent analysis we illustrate a simple, more systematic boundary correction 
procedure in which only the information at the outflow boundary is required. 

The rational for outflow boundary correction is based on the observation that, 
despite the apparently large unsteady region which extends beyond the computa- 
tional domain, the physical source of sound, associated with specific events such as 
the vortex generation and shedding process in the present case or the vortex pair- 
ing in Mitchell et al. (1995a), is captured within the domain. Downstream of the 
vortex-shedding region, the eddies are convected passively and are thus acoustically 
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silent. The boundary correction formula is derived most easily by considering com- 
pact quadrupoles evaluated in two domains of volume integration, whose boundaries 
coincide except at the downstream exits. The exit boundaries, normal to the flow 
direction, are set to be Ayi apart. If Qf : and Q~ are used to denote the first 
time derivative of the Lighthill stress integrals evaluated in the larger and smellier 
domains, respectively, one can deduce that 


Q&t) ~ QiM+.Eii (* - 5) , 

(9) 

53 + Eij (t + -) ■ 

(10) 


The common term in the two equations, Qij(t), represents the physical noise source 
situated within both domains, and function Eij represents the boundary error 
caused by eddies escaping from the respective exit boundary. Since the separa- 
tion between the two boundaries is very small, an eddy can be considered “frozen” 
as it traverses the distance At/i, and hence the error terms in (9) and (10) differ 
only by a small phase difference r = Ayi/17 c , where U c is the local eddy convection 
velocity. Phase-shifting (9) and (10) by r/2 and — r/2, respectively, and subtracting 
the latter equation from the former yield 
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after a Taylor series expansion to 0(r 2 ). The first term on the right-hand side 
of (12) is the algebraic average of the quadrupole sources calculated from the two 
integration domains, whereas the second term clearly represents the desired cor- 
rection at the exit boundary. It is interesting to note that the correction term 
involves only the first time-derivative, in contrast to the second time-derivative in 
the original source terms. In the numerical implementation, Ayj should take the 
smallest possible value (one mesh spacing) to ensure the validity of the frozen-eddy 
assumption. 

Physical insight can be gained by substituting (8) into the right hand side of (12) 
and taking the limit as At/i — > 0, noting that r = Ayi/U c . This leads to 


Quit) = ^2 J v T '> (y» *) d *y +§ij s u * T n (y> o (13) 

which shows clearly that the boundary correction to the quadrupole source calcu- 
lated in a truncated domain Vo is equivalent to the time derivative of the Lighthill 
stress fluxes across the exit boundary So, carried by the convecting eddies at ve- 
locity U c . In other words, the net contribution from the missing acoustic source 
functions outside the integration domain can be approximated by a flux term eval- 
uated on the exit surface. In the above derivation, the eddy convective velocity U c 
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is assumed to be parallel to the j/i axis, and So is a plane surface perpendicular to 
the flow direction. If So is a curved surface, the projected area normal to j/i is used. 
The above formulation can be generalized to the case of noncompact sources. In 
the far-field, the quadrupole contribution to p — 1, denoted as p q , takes the form, 

p ’ (x ’ t] = inr { w X. (JV ‘ ~ Mr) ' ** 

Jr ^J s ^ UcTii{ < 14 ) 

In aeroacoustic problems whose solutions are known to be time periodic, it is 
customary to conduct calculations in the frequency domain (e.g., Mankbadi et al. 
1994; Mitchell et al. 1995a). The corresponding quadrupole noise with a boundary 
correction can be easily obtained by taking the Fourier time transform of (14), 

p 9 (x,u;) = ju> 2 ^ (y,w) e iMur d 3 y 

+iu J s ^UcTij (y, u) c iM “ r <Py | . (15) 

Again, the corrective surface term accounting for the effect of escaping eddies is 
simple and readily computable along with the volume integral. The derivation, on 
the other hand, would be less straightforward if it were carried out in the frequency 
domain where the physical significance of the boundary correction is not as obvious. 

2.8 Model problem 


2.8.1 Near-field simulation 

We consider the unsteady flow field and the sound generated by a NACA0012 
airfoil placed in a 2-D uniform flow at chord Reynolds number Re = 10 4 and two 
angles of attack: a = 5° and a = 8°. In the spirit of Lighthill’s analogy, the acous- 
tic source functions can be determined from am incompressible flow approximation, 
given that the compressibility effect is of 0(M 2 ). A finite-difference code developed 
by Choi (1993) is used to solve the incompressible Navier-Stokes and continuity 
equations numerically in a generalized coordinate system. Second-order central dif- 
ference is used for spatial discretization on a staggered grid. The time advancement 
is of the fractional-step type, in combination with the Crank-Nicolson method for 
viscous terms and the third order Range-Kutta method for convective terms. The 
Poisson equation for pressure is solved using a multi-grid iterative procedure. 

Computations are carried out on a C-type mesh configuration with a total of 896 
by 104 mesh cells. The simulations are run with a time step Af = 2.3 x 10 -3 . No- 
slip velocity conditions are imposed on the airfoil surface. Along the C-shaped outer 
boundary, approximately three chord lengths away from the airfoil, the velocities 
are fixed at the free-stream values, it; = (1,0). At the downstream boundary the 


262 


M. Wang 



x\ 

FIGURE 1. Contours of negative vorticity — u for 2-D laminar flows past a 
NACA0012 airfoil at Re = 10 4 . (a) o = 5°, t — 24.5 (contour levels from —308.0 to 
804.0, increment 8.0); (b) a = 8°, t = 26.8 (contour levels from —335.0 to 1175.0, 
increment 10.0). 

convective outflow condition (Pauley et al. 1987) is applied to allow the vortical 
disturbances in the wake to leave the computational domain smoothly. 

Simulations start with uniform velocity Ui(t — 0) = (1, 0) everywhere. During the 
initial transient period, a starting vortex is shed at the trailing edge and boundary 
layers develop on the upper and lower surfaces. The upper surface boundary layer 
soon separates. It interacts with the lower boundary layer near the trailing edge 
to develop a periodic vortex shedding pattern for a = 5°, as depicted in Fig. la in 
terms of contours of negative vorticity —to at t = 24.5. The same vorticity contours 
for a = 8° at t = 26.8 are plotted in Fig. lb. In this case vortex shedding is initiated 
by the instability of the separated shear-layer near the mid-chord on the suction 
side. The calculated unsteady lift and drag coefficients exhibit aperiodic, perhaps 
chaotic oscillations with time even after an extended time lapse (~ 30 chord flow- 
through times). Similar behavior has been observed and analyzed in the context of 
nonlinear dynamics by Pulliam (1989). 

2.3.2 Acoustic calculation 

The simulated flow-field around the airfoil presented above is two-dimensional, 
implying constant properties in an infinite span. The acoustic formulation in Sec- 
tions 2.1 and 2.2, on the other hand, represents three-dimensional solutions to a 
forced, linear wave equation, and we are interested in the acoustic waves emitted 
from unit span based on the above formulation. Alternatively, one could consider a 
strictly 2-D problem by employing a 2-D version of the acoustic analogy, which can 
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FIGURE 2. The longitudinal quadrupole Qn calculated from three different-sized 

source domains whose downstream boundaries are located at x\ ~ 2.59 ( ), 

2.89 ( ), and 3.21 ( ). The airfoil angle of attack a = 5° (a) Without 

boundary correction; (b) with boundary correction. 


be derived easily by integrating the 3-D formulae along the infinite span (Mitchell et 
al. 1995b; Wang 1993). The results are, however, of less physical relevance since in 
practical situations involving a long span, the near-field inevitably develops three- 
dimensionality, and thus the phase difference between the various radiating elements 
along the span cannot be ignored. 

To illustrate the effect of boundary correction, Figs. 2a and 2b compare the time 
oscillations of the compact longitudinal quadrupole Qn calculated from (13), before 
and after the boundary correction term is added. In the calculations, it is assumed 
that Tij « ViVj, an approximation justified by the relatively large Reynolds number 
and small Mach number. The three curves represent evaluations based on three 
different sized source domains whose downstream boundaries are 20 grid points 
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FIGURE 3. The same as in Fig. 2, except that the airfoil angle of attack a = 8°. 


apart, located at y\ « 2.59, 2.89, and 3.21, respectively. The primitive (uncorrected) 
Lighthill quadrupole, shown in Fig. 2a, is seen to exhibit a strong dependence on the 
downstream boundary location. After applying the corrections, the three curves are 
seen to converge as shown in Fig. 2b, indicating the physical noise source caused by 
vortex generation near the trailing edge, which is well captured within all the three 
integration domains. In estimating the boundary Reynolds stress fluxes, a constant 
convective velocity U c = 0.90 is used for all the three surfaces. The corrective effects 
are found to be rather insensitive to the value of U c , although a slight improvement 
has been observed by fine-tuning the value of U c based on the local convective 
velocity on each surface. The other two quadrupole components, Qu and Q 22 , 
are computed and compared in the same manner, and equally drastic reductions 
in boundary errors are obtained. The residual boundary error for Q 22 is, however, 
larger than those in Fig. 2 due to its larger pre-correction error magnitude. 

An example of boundary corrections applied to an aperiodic source is given in 
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Figs. 3a and 3b. The quadrupole Qu is obtained based on the near-field simulation 
data for the case of a =8°. Again, three exit boundaries passing the same yj 
stations as in Fig. 2 are used to illustrate the dominant impact of boundary errors 
in the source calculation, as shown in Fig. 3a. The corrected Qn source terms in 
Fig. 3b, as well as the other two quadrupole components not shown in the figure, 
compare well despite the non-harmonic nature of the signals. This is important 
since the main advantage of solving aeroacoustic problems in the time domain lies 
in its ability to treat arbitrary, non-periodic signals. In comparison with the a = 5° 
case (cf. Fig. 2b), Fig. 3b indicates higher amplitude, lower frequency acoustic 
phenomena corresponding to the stronger vortices shed at a slower rate in the near- 
field. The eddy convective velocity U c is found to be best approximated by 1.0. 

In the event that the source integration domain is not small in comparison to 
the dominant acoustic wavelength, retarded-time variations in the source region 
become significant. The boundary corrections are applied directly to the far-field 
density according to (14). To examine the efficacy of this more general approach, 
an efficient integration-interpolation method has been developed for the evaluation 
of the surface and volume integrals in (5) and (14). The scheme treats each near- 
field computational cell or boundary element as an individual acoustic source. At 
each simulation time step t n , the future time a n = t„ + M |x — y|, at which the 
emitted acoustic signal reaches the far-field position x, is calculated. The scheme 
then locates in the discretized far-field time series the point straddled by a n -i and 
cr n , and interpolates linearly on the integrands to find their contributions to the 
far-field density at that time. The total density history at the observation point is 
obtained by summing up contributions from all the source elements and simulation 
time steps. This procedure, which is second order accurate, uses the same time step 
as for the near-field simulation and requires minimal extra computer memory. 

Figures 4a and 4b contrast the far-field acoustic pressure signals due to quadrupole 
radiation for the case of a = 5°, evaluated in the same integration domains as 
in Fig. 2, before and after the boundary corrections. The fluctuating pressure 
(p q = 7 p q , renormalized by the mean free-stream value) is calculated directly from 
(14) with full account of the retarded time, at a given far-field position |x| = 50, 
0 = 30° , where 6 defines the angle measured counter-clockwise from the downstream 
x\ axis. The free-stream Mach number M = 0.2. The three p q curves correspond- 
ing to the three different exit boundaries compare rather well after corrections are 
applied. The discrepancy, in the form of higher frequency oscillations, is mainly due 
to the T22 source component whose boundary effect is more difficult to eradicate. 

The retarded-time effect is demonstrated in Figs. 5a and 5b for M = 0.1 and 
0.2, respectively, by comparing the computational results with and without using 
the compact source approximation. The two figures again plot the quadrupole 
contribution to the acoustic pressure at |x| = 50 and 0 = 30°, for the case of 
a = 5°. The solid lines are obtained using (14), while the dashed lines are based 
on the compact source formulae (6) and (13). Both have been subjected to the 
appropriate boundary corrections. The integration domain employed is the smallest 
among the three used for Fig. 4. As expected, the noncompact formulation produces 
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FIGURE 4. The acoustic pressure at |x| = 50, 6 = 30° due to quadrupole radiation, 
calculated from (14) using the same three integration domains as in Figs. 2 and 3. 
The airfoil is flying at M = 0.2 and a = 5°. (a) Without boundary correction; (b) 
with boundary correction. 

stronger radiation because of less cancellation among signals from different source 
elements. The retarded-time effect is relatively small at M = 0.1 and becomes 
somewhat significant at M = 0.2. Furthermore, one observes that there exists little 
phase difference between the compact and noncompact solutions, suggesting that 
the acoustic source is centered near the origin of the far-field coordinate x, which 
in the present calculation is defined at the trailing edge. 

It should be pointed out that in the examples described in Figs. 4 and 5, the 
active source (vortex-shedding region) is only slightly noncompact. On the other 
hand, the source integration domains, containing 6 to 10 eddies depending on the 
exit boundary selection, are fairly noncompact, particularly for the case of M = 0.2. 
They thus provide valid tests for the general boundary treatment method (14) as 
well as the numerical integration scheme with retarded time variations. The good 
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FIGURE 5. The acoustic pressure at |x| = 50, 0 — 30° due to quadrupole radiation, 

calculated based on noncompact ( ) and compact ( ) source formulations. 

The airfoil angle of attack a = 5°. (a) M = 0.1; (b) M = 0.2. 


agreement with compact source solution at low Mach number (M = 0.1 and below) 
shows that the integration procedure is capable of providing adequate cancellations 
among signals from the passive convection region. 

A comment is in order concerning the selection of the integration boundary So- 
In principle, So can be placed anywhere so long as the active noise source is enclosed 
within the source integration domain. From a computational standpoint, however, 
the velocity field in the vicinity of the outflow boundary for the Navier-Stokes 
simulation is somewhat distorted due to the application of the convective boundary 
condition. The non-physical eddy distortion may serve as another source of spurious 
noise if it is included in the calculation. For this reason, we always choose smaller 
source integration domains than the actual flow simulation domain by placing So 
at least 15 grid points upstream from the computational outflow boundary. 

Finally, the acoustic dipole sources due to the unsteady compressive stress exerted 
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FIGURE 6. Acoustic dipoles on the airfoil surface, calculated from (8) for (a) 
a = 5°; and (b) a = 8°. D\ (drag dipole); I >2 (lift dipole). 


by the airfoil surface are evaluated and compared with the volume quadrupoles. Un- 
der the condition that the airfoil chord is small relative to the acoustic wavelength, 
the compact solution form (6)-(8) applies. The calculated dipole components are 
depicted in Figs. 6a and 6b for the cases of a = 5° and a = 8°, respectively. Both 
the pressure and viscous stress are included in the calculations of pij. As expected, 
the viscous contribution to the lift dipole D 2 is negligibly small. Its effect on the 
drag dipole D \ , on the other hand, reaches approximately 27% in magnitude with 
opposite phase, relative to the pressure contribution, for the case of a = 5°. Overall, 
the lift dipole is much stronger than the drag dipole, and both amplify with increas- 
ing angle of attack. Relative to the quadrupole components calculated earlier, £>2 
is of larger or comparable magnitude and thus dominates the far-field radiation, 
given that its coefficient in (6) is 0(M _1 ) larger than that of the quadrupole terms. 
Along the horizontal (ri) axis, however, the quadrupole sources may play a limited 
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role because the magnitude of D\ is also small. 

3. Conclusions and future directions 

In acoustic analogy based calculations of aerodynamic sound, it is important to 
distinguish the active source region, where production of the unsteady Reynolds 
stress takes place, from passive regions characterized by the convective motion of 
eddies. The former should always be enclosed within the domain of source integra- 
tion, whereas the latter may be truncated provided that an adequate account of the 
effect of eddies crossing the permeable integration boundary is provided. 

This report illustrates the necessity for, and means of, boundary corrections 
through a paradigm problem of airfoil vortex shedding. It demonstrates that the 
spurious noise generated by the exit boundary in the airfoil wake is due to the 
time variation of the unsteady momentum fluxes across the boundary, carried by 
the escaping eddies. For a class of problems where eddies (organized structures) 
leave (enter) the source integration domain at a nearly constant speed, the spu- 
rious boundary noise can be eliminated or reduced drastically by the correction 
terms in (12)-(15). This approach allows a quantitative evaluation of the radiated 
quadrupole noise, which has rarely been done in the past. 

In the case of vortex-shedding noise from a small airfoil, studied above, compu- 
tational results suggest that the volume quadrupole noise is small in the low Mach 
number limit in comparison to the lift and drag dipole noise emitted from the airfoil 
surface. The techniques developed in this study are equally applicable to other flow 
configurations such as jets and mixing layers, where the volume acoustic sources 
provide the dominant contribution to the far-field noise. 

A crucial issue to be addressed next is the aeroacoustic scattering by the air- 
foil trailing-edge, the major source of broadband noise according to experimental 
measurements (Brooks &: Hodgson 1981; Blake & Gershfeld 1988). Both theory 
(Crighton & Leppington 1971; Howe 1978) and experimented results indicate a 
non-multipole character of the sharp-edge noise. In the limiting case of scatter- 
ing by a half-plane, the far-field intensity exhibits a M 5 dependence and sin 2 (0/2) 
directivity. A computational prediction should capture these fundamental charac- 
teristics. While the Curie integral is formally exact, its usefulness as a predictive 
tool depends upon the precise knowledge about the volume and surface source 
terms. If the near-field source functions are approximated by the incompressible 
Navier-Stokes solutions, as in the present case, the Curie-integral prediction is only 
accurate for compact surfaces relative to the dominant acoustic wavelength, since 
the exact boundary conditions on the surface are not satisfied by the acoustic field. 
For acoustically noncompact airfoils where edge scattering noise is generated, the 
compressible (acoustic) contribution to the surface integral is significant, and proper 
account of the acoustic-surface interaction must hence be taken. 

In order to account for the surface reflection and edge diffraction, the appropriate 
Green’s functions, whose normal derivatives vanish on the solid surface, must be 
employed in an integral solution to the Lighthill equation. Solution development 
may be carried out in either the time domain or the frequency domain. In either 
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case, it is likely that the Green’s function has to be determined numerically in 
the discrete sense, given the relatively complex geometry. An alternative to the 
Green’s function approach would be to solve the Lighthill equation numerically, 
using the turbulence data generated by the near-field Navier-Stokes simulation as 
a discretized forcing function. These different approaches will be explored to find 
the most effective solution technique. In the meantime, we will work closely with 
the LES group at CTR in an effort to establish reliable turbulence statistics in the 
near field to be used as source functions for the acoustic analogy based far-field 
calculations. 
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